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ABSTRACT 

Observations of nearby molecular clouds detect “shells”, which are likely caused by winds from young 
main sequence stars. However, the progenitors of these observed features are not well characterized 
and the mass-loss rates inferred from the gas kinematics are several orders of magnitude greater than 
those predicted by atomic line-driven stellar wind models. We use magnetohydro dynamic simulations 
to model winds launching within turbulent molecular clouds and explore the impact of wind properties 
on cloud morphology and turbulence. We find that winds do not produce clear features in turbulent 
statistics such as the Eourier spectra of density and momentum but do impact the Eourier velocity 
spectrum. The density and velocity distribution functions, especially as probed by CO spectral lines, 
strongly indicate the presence and influence of winds. We show that stellar mass-loss rates for indi¬ 
vidual stars must be rh^ > I0“^MQyr“^, similar to those estimated from observations, to reproduce 
shell properties. Consequently, we conclude that B and A-type main sequence stars have mass-loss 
rates several orders of magnitude larger that those predicted by models or that young stars are more 
variable than expected due to magnetic activity or accretion. 

Subject headings: stars: formation, ISM: bubbles, ISM: clouds, stars:winds, outflows, ISM: jets and 
outflows, turbulence, magnetohydrodynamics (MHD) 


1. INTRODUCTION 

Turbulence within molecular clouds is thought to be an 
essential component of the star formation process, influ¬ 
encing cloud lifetimes, reducing star formation efhcien- 
cies and even setting the distribution of ste llar masses 
(|McKee & Ostriker||2007t jOffner et al.||20I4| and refer- 

molecular cloud 


ences therein). However, the origin of 
turbulence and the means by which it is regenerated re¬ 
main debated. One possibility is that young stars them¬ 
selves inject sufficient energy to drive turbulent motions 
through protostellar out flows or radiatively driven winds, 
at least on parsec scales (Erank et al. pMil and references 
therein). Signatures of such activity exist within molec¬ 
ular clouds, however, it is difficult to quantify the true 
impact of these sources on the cloud motions. 

Large-scale surveys of the Galactic plane, such as the 
GLIMPSE survey performed with the Spitzer Space Tele¬ 
scope, have discovered thousands of shell-like features in 
which molecular gas has been swept up and evacuated 
from large regions. Many of the shells are directly asso¬ 
ciated with high-mass 0-type stars, which produce ioniz- 
ing photons as we l l as en ergetic radiatively driven winds. 
Deharveng et al. (2010) investigate a sub-sample of the 
Spitzer Galactic plane data and And that ^ 15% of the 
shells do not enclose ionized gas, which implies they may 
be driven by A or B-type stars. 

Recently, high-resolution spectral observations of the 
nearby P erseus molecular cloud have identified a number 
of shells (|Arce et al.|2011| henceforth All). The authors 
estimate that the momentum and energy in such features 
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could account for approximately half of the total cloud 
turbulent energy, much more than could be contributed 
by collimated outflows from protostars. The larger shells 
are expanding with velocities in excess of the typical tur¬ 
bulent cloud velocity and create clear kinematic signa¬ 
tures in the gas. The shell regions influence the cloud ve¬ 
locity distribution significantly such that it deviates from 
the turbulent velocity distribut ion predicted by cloud 
simulations excluding feedback (Beaumont et al. 2013). 
Thus, it seems possible that the statistical properties of 
the cloud turbulence, which cascades from the largest 
scales, are appreciably different for turbulence influenced 
by internal stellar sources. Eor example, there is some 
evidence that turbulence metrics such as the momentum 
and energy Eourier spectra are di fferent for protostel¬ 
lar outflow driven turbulence (e.g.. Swift & Welch||2008 
Carroll et al.||2QIQ ). 


Within lower mass star forming molecular clouds like 
Perseus, 0-type stars are not present. Instead, some of 
the shells appear to originate from less massive, but more 
numerous, B-type stars. Several such candidate progeni¬ 
tors are observed in and around the Perseus cloud (All). 
However, a number of the progenitors do not appear to 
have either associa ted O or B type stars a nd the shell 
origin is unknown. Nakamura et al. (2012) And parsec- 
scale bubbles surrounding two small star clusters in the 
Orion region. They suggest that the bubbles are driven 
by the combined winds from the most massive stars in 
these clusters, which are B and A-type stars. This could 
indicate that winds from multiple lower mass stars in a 
cluster may also inject significant energy to create such 
features. The B-type main sequence star p Ophiuchi, 
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which is located near the Ophiu chus star forming; reg ion, 
also has a sizable shell of 1.4pc (jWilking et al.||2008). 

If A and B-type stars generate some of these 
shells through radiatively driven winds, this presents 
a quandary for stellar atmospheric models. Mass-loss 
rates estimated from the shell masses, gas velocities, 
and Perseus cloud properties suggest winds with ^ 
10-^ - IQ-^M© yr “h This is more than three orders 
of magn itude higher than those predicted by theoretical 
models ( Smith|2Q14 see their Figure 3). In fact, around 
the transition mass between O and B-type stars, radia¬ 
tive models suggest that the stellar radiation field may 
be too weak to drive mass-loss through optically thick 
spectral lines; this is known as the “weak wind prob¬ 
lem”. Even without the complication of the weak wind 
problem, mass-loss rates for lower mass stars extrapo¬ 
lated from higher mass values would rem ain orders of 
magnitude too small (e.g. Vink et al.||2QQl ). 


To date, little numerical work has been devoted to 
modelling winds from intermediate mass stars. Prior 
work has focussed exclusively on winds from higher mass 
0 -ty pe stars and their impact on the gas in young clns- 


ters (Dale & Bonnell 2008 Rogers & Pittard 2013). In 
these cases ionizing radiation plays a key role in cre- 


(IGendelev & Krumholz 20121 Dale et al. 

1 1 , n — 

20131 Rogers 

& Fittard 

im 

IGeen et al.l 20151. Only 

Dale ^ Bon- 

nell ( 2 UU 8 

) and 

Dale et al. (|2013|) have considered pure 


model mass-loss from stars with m > lOM© and do not 
include magnetic fields. 

In this paper, we use hydromagnetic simulations to ex¬ 
plore the range of mass-loss rates needed to reproduce the 
properties of observed shells in local star-forming clouds. 
We apply a variety of statistics to assess the degree to 
which the winds impact gas motions and assess whether 
winds impart detectable signatures. We construct syn¬ 
thetic observations of the simulations in CO( 1 - 0 ) to as¬ 
sess the accuracy of observational metrics in identifying 
and estimating properties of observed shells. We describe 
our simulations and models in ^ We present results, in¬ 
cluding turbulent statistics and synthetic observations in 
^ We discuss the results and evaluate alternative ex¬ 
planations in Qand summarize our conclusions in ^ 

2. SIMULATION PARAMETERS 
2 . 1 . Numerical Methods 

For our study we use the magneto-hydro dyna mics 


adaptive mesh refinement (AMR) code ORION (e.g., My- 
ers et al.|2014|). ORION uses the chombo toolset to solve 
partial ditterential questions on block-structured AMR 
meshes. In the simulations we employ the constrained- 
transport m agneto-hydrodynamic (CTMHD) module 
et al.| 2012 ), which is base d on the finite volume schem e 
implemented in PLUTO (Mignone et al.| 2Q Q7[ 2012 ). 
ORION also has Lagrangian sink particles (|Krumholz 
et al.||2004|), which are inserted when collapsing gas ex¬ 
ceeds tHeJ cans criterion on the maximum AMR level. 
Extensions of the sink particle methodology allow ORION 
to compute protostellar luminosities, launch protostellar 
bip olar outflows, and handle the accretion of magnetized 
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Fig. 1. — Distribution of observed (black solid line) and simulated 
(colored lines) mass-loss rates. WI and W2 indicate the mass-loss 
distributions for the models shown in Table □ 


functionality: a sub-grid model that allows ORION to fol¬ 
low the launching of isotropic stell ar w inds. We present 

ORION also has 


the wind methodology details in §2.3 
modules to compute self-gravit y and radiative transfer 


via grey flux-limited diffusion (Miniati & Colella 2007 
Krumholz et al.j|2007 ), tiowever, we do not include these 
in the present simulations. 

We treat the gas as an isothermal ideal gas such that 
its thermal energy is e = pksT/ [fi{j — 1)], where p is 
the gas density, is the Boltzmann constant, T is the 
gas temperature, p is the mean particle mass, and 7 is 
the ratio of specific heats. We adopt 7 = 1.0001 and 
p = 3.9 X 10“^^g, which is the typical local value for cold 
molecular gas comprised of H 2 and He. 


2 . 2 . Global Parameters 

The simulations model a piece of a molecular cloud 
using a box with side L = 5 pc and M = 3762 M©. 
The initial cloud temperature and Mach number are 
T = 10 K and MIsd = 10.5, respectively, which places 
the simulated cloud on the line-width size relatio n, ctid = 


McKee & ©striker 


2007). To ob- 


gas ([Offner et al. 2009 Cunningham et al. 2011 Lee et al 

20141). In this paper we implement a new sink particle 


0.72Rp^^ kms ^ (e.g., 
tain initial density and velocity conditions, we “drive' 
the simulation gas for approximately two mach crossing 
times without gravity by adding ra ndom large-scal e per¬ 
turbations to the velocity field (e.g. [Mac Low|l999| ). The 
perturbation wave numbers are between 1 and 2 , which 
corresponds to a linear size oi L/2 — L. 

The simulations start with an initially uniform verti¬ 
cal magnetic field, which becomes disordered over 
the turbulent driving period. The simulations have 
plasma parameters (ratio of thermal pressure to mag¬ 
netic pressure), P = S7TpQcllB‘^ = 0.02 — 0 . 6 , where 
Cs is the sound speed and po = 2.04 x 10“^^g cm“^ is 
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TABLE 1 

Model Properties^ 


Model 

B{nG) 

/9 

ti (Aross) 

hun(Myr) 

Mtot(lO-6M0yr-l)b 

'i^max (kms 

(AU) 

Driving 

W1_T1_L0 

13.5 

0.1 

1.6 

0.15 

41.7 

200 

4x10^ 

N 

W1_T1 

13.5 

0.1 

1.6 

0.1 

41.7 

200 

10^ 

N 

W0_T2 

13.5 

0.1 

2.0 

0.2 

0 

200 

10^ 

N 

W1_T2 

13.5 

0.1 

2.0 

0.2 

41.7 

200 

103 

N 

W2_T2 

13.5 

0.1 

2.0 

0.2 

4.5 

200 

103 

Y 

W2_T3^ 

5.6 

0.6 

2.0 

0.1 

4.5 

200 

10® 

Y 

W2_T4‘^ 

30.1 

0.02 

2.0 

0.1 

4.5 

200 

10^ 

N 

W3_T2d 

13.5 

0.1 

2.0 

0.2 

1 X 10-3 

103 

103 

Y 

W4_T2 

13.5 

0.1 

2.0 

0.15 

1 X 10-2 

103 

103 

Y 


^ Model name, initial mean magnetic field, ratio of thermal to magnetic pressure, the initial start time in crossing 
times, the run time, the total stellar mass-loss rate, the maximum wind velocity, minimum grid resolution, and 
whether external driving is continued simultaneously with winds or whether winds alone inject energy. All models 
have L = 5pc, M = 3762Mq, Ti = lOK and = 5. All the calculations are first evolved with external driving but 
without sources to set up the initial turbulent conditions at t = 0. 

^ The estimated mass-loss rate from all stellar winds in Perseus is 9.49 x 10“®MQyr“^ (Arce et al. 2011). 

Same wind model as W2_T2 but a different initial magnetic field. 

^ Same initial turbulence as run W2_T2, but containing only a single isolated star with the given mass-loss rate. 


the mean density. The fiducial Alfven Mach number is 
= \/127r/)oc^iD/^o = 2.3. The simulations follow 
four different wind launching models (W1-W4) and four 
different initial turbulent patterns (T1-T4), which share 
the same Mach number but vary spatially and/or with 
bulk magnetic field strength at the initial time. The sim¬ 
ulation parameters are summarized in Table 
In W1 and W 2 , we randomly place five stellar sources 
with different mass-loss rates. W1 has mass-loss rates 
ranging from 5.7 x 10“^ — 2.0 x 10“^MQyr“^ and W2 
has mass-loss rates ran^^g from 2.6 x 10 “^ — 2.0 x 
lO“^M 0 yr“^ (see Figure^^. This number of sources is 
equivalent to the source number per volume identified 
in Perseus. The sources are meant to represent main se¬ 
quence stars that have either formed in an earlier genera¬ 
tion of stars or are older stars that formed elsewhere and 
wandered into the region. Either of these scenarios pro¬ 
vides a plausible explanation for the observed sources in 
Perseus. For example, a 20 pc size cloud like Perseus has 
a dynamical time of tdyn — LIcjid ^ 9 Myr, which is suf¬ 
ficient time to form stars of 3-5 solar masses, whic h would 
reach the main s equence in 0.6 — 3 Myr (e.g. Waller- 


stein et al. 1997). Alternatively, the observed sources 
may have wandered into Perseus from a nearby region. 
B-type stars remain on the main sequence for approxi¬ 
mately 1 — 5 X 10 ^ years. If a star ha s a typical velocity 


Foster et al. 


2015), then it can travel 


of 1 kms“^ (e.g., 
about 1 parsec in 1 IVIyr, or a tew hundred parsec before 
the end of its main sequence lifetime. This distance is 
sufficient such that some of the observed sources could 
have formed hundreds of millions of years earlier in a 
nearby cloud, whose gas is now dispersed. 

In two additional simulations, W3 and W4, we exam¬ 
ine the wind evolution for a single star with a prescribed 
mass-loss rate of 10“^ or 10“^ M 0 yr“^. These values are 
the mass-loss rates predicted for earl y O and late- t ype B 
stars, respectively, according to the Vink et al. (2001) 
model. Following the evolution of an isolated source 
avoids interactions between winds, which occurs in the 
other simulations and may boost the efficacy of gas en¬ 
trainment and energy injection. 

Finally, we also compute a purely turbulent model with 
no wind sources, WO. This model has the same initial 


turbulence as model T2, but the turbulence is simply 
allowed to decay without any energy input. 

During the wind phase, we add additional refinement 
to the simulation s according to the T ruelove criterion 
with Nj = 0.125 ( Truelove et al.||1997 ). The application 
of a density gradient criterion, Ap/p < 10 , ensures that 
wind shells are resolved by two additional levels of refine¬ 
ment for a maximum resolution of Ax 2 = 10^ AU. We 
also rerun one simulation with no additional refinement 
to check the sensitivity of the results to the grid resolu¬ 
tion (W1_T1_L0). In some of the models, once the stars 
are inserted and wind launching begins, external turbu¬ 
lent driving ceases and winds provide the only energy 
injection. Whether driving is continued after the initial¬ 
ization phase is indicated in Table In other cases, the 
models have both winds and external driving. 


2.3. Wind Model 

The three key properties of the wind model are the 
mass-loss rate, the wind velocity, and the wind tempera¬ 
ture. Stellar mass is a factor only insofar as the mass-loss 
rate depends on it, otherwise, mass has no direct impact 
on the calculation since it has no dynamical impact with¬ 
out gravity. Star masses are randomly drawn such that 
their mass-loss rates are comparable to those inferred 
from the observed sources, i.e. 3 M 0 < < 15M0. 

The most fundamental parameter of the simulation is 
the total mass-loss of the stars, Mtot- A random deter¬ 
mination of 5 sources is inherently very stochastic, and a 
single object may dominate the total mass-loss budget. 


2.3.1. Mass-loss Rate 

The mass-loss rate for a given star can be related to 
fundamental param eters via an analyt ic fit to observed 
data as presented in |Vink et al.| (|2001|) . The fit depends 
on five stellar pararneters: the stellar luminosity (T^), 
stellar mass (m*), effective temperature (Teff), metallic- 
ity (Z), and L’rat = '^inf/'^esc? which is the ratio of the 
terminal velocity to the escape speed at the surface. For 
Teff < 2.25 X lO^K: 


log(m) = — 6.688 + 2.210 log 



- 1.339 log 
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-I.OOIl„s(||)+ 1.07 log 

+0.85 log{Z). (1) 

For Teff > 2.25 x lO^K: 


log(m) 


-0.697+ 2.104 log (+^)-1.313 log ( 5 ^) 
-1.226 log (II)+0.933 log 


-10.92 


log 


TeS 


4.0 X lO^K 


2 

+ 0.85X log(Z). 


( 2 ) 


Prom Lamers et al. (1995): 


ti,^t = 0.7 Teff < 1.25 X 10"^ K 

Vrat = 1-3 1.25 X 10^ < Teff < 2.1 X 10^ K 

Vrat = 2.6 Teff > 2.1 X 10^ K. 


( 3 ) 

( 4 ) 

( 5 ) 


Note the values of the effective temperature for the tran¬ 
sitions in Vj-at are somewhat approximate, and different 
authors adopt slightly different values. 

Although the stellar masses, and thus the mass-loss 
rates, are time dependent, the change is small over the 
simulation evolution, so the mass-loss rate can effectively 
be treated as a constant. Figure shows the distribu¬ 
tion of mass-loss rates for the observations and simulated 
sources. 


2.3.2. Wind Velocity 

The terminal wind velocity for OB line dri ven winds is 
a few hundred to a thousand kms“^. From ivink et al.l 

( |20M] ): 

log(u) = 1.23 — 0.3 log(L) -f 0.55 log(m) -f 0.64 log(Teff) 
-1-0.13 log(Z). (6) 

In the simulations, we limit the wind launching speed 
to some maximum velocity, Umax (see Table , in order 
to avoid inefficiently small time steps. Lower mass stars 
have smaller mass-loss rates and higher velocity winds. 


evidence of dust emission in some of the bubble interiors, 
which would require temperatures T < 1100 K, the dust 
sublimation temperature. 

The thermal pressure of the wind has little impact dy¬ 
namically. For example, if the thermal pressure of a gas 
with a average temperature, T, and density, p, is given 
by Pth = where Cg = y/ksT /then the fidu¬ 
cial simulation parameters yield an average molecular gas 
pressure: 


P/kB = 5224 f- V. -?i ( i Kcm“^ 

\2.04 X 10“^^ gcm“^y \10 Ky 

( 7 ) 

The pressure inside a wind bubble is the sum of the ther¬ 
mal pressure and gas ram pressure, where both depend 
upon the total mass-loss: 


-^wind -Pram 4“ Pbh 

TTlwt 2 2 

= -7; - C 

V ^ V 


( 8 ) 

( 9 ) 


where V is the bubble volume. In a turbulent medium 
the evacuated region will not be exactly spherical, but 
we can still approximate V (x R^. For a lO^K wind 
with terminal velocity v = 200 kms“^ the ratio of the 
ther mal pressure to th e ram pressure is 1.5 x 10“^ (see 
also Geen et ar]|2015 ). 


2.4. Relation between Mass Loss and Shell Evolution 

We can use a simple calculation to relate the stellar 
mass-loss rate to the shell expansion. At all times the 
mass in the shell will be dominated by the swept up ISM 
mass not the integrated wind mass from the star. Thus, 
we can define the total shell mass as 

Ms = -R^Po, ( 10 ) 

where po is the average cloud density (pre-wind), Rs is 
the shell radius, and we assume the shell is roughly spher¬ 
ical. The expansion momentum of the shell is given by 

MigVg — ^ (^^) 


2.3.3. Wind Temperature 

We assume that the wind is ionized, and thus set the 
launched gas temperature to 10^ K and compute its ther¬ 
mal energy using p = 2.13 x 10“^^g representing a mix 
of atomic H and He. We do not include cooling explic¬ 
itly, such that hot wind gas can only cool by mixing 
with ambient cold material. Although lower mass stars 
like the Sun have coronae that are cooler, their winds 
are heated to much higher temperatures by convection 
( Owocki|2QQQ ). Conversely, the winds of high-mass stars, 
which do not have a convective zone, are similar to that 
of their corona (10^ — 10^ K). 

Inspection of the Perseus CO maps indicates that there 
is no CO inside most (but not all) of the bubbles. This 
suggests that the molecular gas has been dissociated and 
the temperatures inside the cavities are above ~ 10^ K. 
H+ is detected in one or two of the bubbles, namely 
the ones that appear to be driven by higher-mass stars. 
There is likely to be atomic hydrogen in some of the 
bubbles, but HI maps with sufficiently good resolution 
to confirm this directly do not exist. However, there is 


where Vs is the shell expansion velocity, rh^ is the mass- 
loss rate, is the wind velocity and t^ is the age of the 
shell (length of time the wind is driven by the star). This 
time, can be related to the size and expansion velocity 
of the shell. At early times, the expansion follows the self¬ 
similar Taylor-Sedov blast wave solut ion (radiative losses 

are negligible): Rg oc tw^ (Shu 1992). This corresponds 
to a shell velocity Vg = 2/b[Rg/tw)^ which gives: 


2R 

MsVs = 


Then the mass-loss rate is: 

l^ViRlpo 


rUn 


?>Vu 


( 12 ) 


(13) 


where the mass of the shell has been replaced using Equa¬ 
tion [Tol The mass-loss rate can be normalized in terms 
of fiducial values: 


rhw = 5.24 X 10 


-7 




1 kms 


200 kms 


-1 


Rs 
1 pc 
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Fig. 2.— Volume rendering of the gas velocity, where gray roughly traces the ambient cloud and red highlights the high-velocity wind 
regions. Top: Run W2_T2 (/3 = 0.1) at times 0 Myr and 0.1 Myr (right). Top-Middle: Run W1_T1 at times 0 Myr and 0.1 Myr (right). 
Bottom-Middle: Run W2_T3 (/3 = 0.6) at 0.1 Myr and Run W4_T2 (M = 10“® Mq yr“^) at 0.15 Myr (right). Bottom: Run W3_T2 
(M = 10“^ MQyr“^) at times 0.1 Myr and 0.2 Myr (right). Runs W2_T2, W3_T2 and W4_T2 have the same gas distribution at t = 0. 
For Runs W3_T2 and W4_T2 the wind source is in the center of the image. 
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_3 )M0yr \ 
\200 cm J 


(14) 


where no is the mean particle number density of the cloud 
^0 = Here, 1km s“^ is a typical measured shell 

expansion velocity and 1 pc is a typical shell radius. 

3. RESULTS 

We analyze the simulation outcomes by computing 
the bulk cloud properties, turbulence statistics, and syn¬ 
thetic observations of the clouds. 

3.1. Bulk Properties 

Figure shows volume rendered images of the gas ve¬ 
locity for the models at different times. The regions 
of high velocity gas are roughly centered on the stellar 
sources and are co-spatial with the hot 10^ — 10^ K gas. 
The shells are more spherical at early times and loose 
some symmetry as they expand into lower density re¬ 
gions. The magnetic field has a confining and stabilizing 
effect, which suppresses shell instabilities that might oth¬ 
erwise occur at t he interaction boundary between dense 
and light flui^ (|Krumholz et al. 2009 Arthur et al. 
2011|). Figure |3| shows the magnetic held vectors overlaid 
on slices of the gas density at four star locations. The 
shells impact the local magnetic field morphology such 
that the field is mo stly parallel to the shell boundaries. 
Ridge et al.| (|2006a|) note a similar configuration for the 
shell near the i(J348 cluster in Perseus, where polariza¬ 
tion measurements indicate the magnetic field aligns with 
the shell boundary. As gas is swept up by the expanding 
shell, the bulk magnetic field, which is coupled to the gas, 
also declines within the evacuated region. Figure also 
illustrates that the shell morphology is a strong function 
of the local gas density and turbulent distribution. 

A close examination reveals that the shell walls are 
not completely smooth. A few tendrils of gas on the 
shell boundaries protrude into the hot, low-density re¬ 
gion (e.g., see the right portion of the shell in the top- 
left panel). These features increase with magnetic field 
strength, leading to noticeable wiggles in the shell sur¬ 
face, and disappear in the weakest field case. We con¬ 
clude they are likely caused by a magnetic kink i nsta¬ 
bility, which also develops in astrophysical jets (e.g.. 


Lee 


et al. 2015). We might expect to see Rayleigh-Taylor 


instability, which frequently occurs at the interfaces be¬ 
tween low-density and high-density fluids, or Kelvin- 
Helmholz instability, which results from she ar between 
two fluid layers (e.g., Ntormousi et al. 2011). However, 
we do not see any clear evidence of either instability. This 
may be because the gas is turbulent or because they are 
suppressed by the magnetic field. 

Runs W3_T2 and W4_T2 follow the evolution of in¬ 
dividual isolated stars with mass-loss rates of 10 “^ and 
10“^ M 0 yr“^, respectively. These winds have such a 
minor sphere of influence that the high-velocity bubble 
region appears only very faintly in the middle and bot¬ 
tom panels of Figure As illustrated in Table each of 
the wind models have very different hot gas filling frac¬ 
tions, which correlate roughly with the total mass-loss 
rate. 

Figure]^ shows the evolution of the total velocity dis¬ 
persion, rms magnetic field strength and effective size 
of the bubble region. Only the highest mass-loss rate 
case (Run W1_T1 ) shows appreciable energy injection 


by the winds, as indicated by the 50% increase in the 
gas velocity dispersion. However, after an initial phase 
of rapid growth, the shells expand more gradually. The 
shell around the source with the highest mass-loss gen¬ 
erally dominates the bulk properties. 

Run W1_T1 is the only run which exhibits a notice¬ 
able change in the rms magnetic field. The magnetic 
field increases because the field is swept up and com¬ 
pressed by the shell expansion. The diffuse gas in the 
bubble interiors have negligible magnetic field. Figure 
shows that the evolution of the runs with and without 
AMR refinement (run W 1 _T 1 and run W1_T1_L0) are 
nearly identical, indicating that these results are fairly 
independent of the shell resolution. 

In the lowest mass-loss cases, as illustrated by runs 
W3_T2 and W4_T2, the bubbles rapidly reach pres¬ 
sure equilibrium with the surrounding gas and expan¬ 
sion stalls. Figure [^indicates the time for radiative bub¬ 
bles to be c ome p ressure-confined as predicted by |Koo| 
& McKeel ([19921, where we modify the prediction to 


account for turbulent and magnetic pressure (see Ap¬ 
pendix®. We find the time at which expansion stalls is 
in goodagreement with the predictions for the mass-loss 
rates. Runs W3_T2 and W4_T2 expand slightly longer 
since their peak wind speed is ^ 600 — 700km s“^ rather 
than 200 km s“^, which gives a slightly longer confine¬ 
ment for their mass-loss rates. In the lowest mass-loss 
case {ttIw = 10 “^ M 0 yr“^) the bubble volume increases 
and decreases in size depending on fluctuations in the 
local turbulence and gas density. Between 0.1 Myr and 
0.2 Myr the shell becomes so small that it is invisible on 
the sc ale of the image in Fi gure ^ 

The Koo & McKee (1992) mood predicts that a source 
with = 10 “^ M 0 yr“^ will become pressure confined 
at 0.22 Myr for our fiducial parameters, which means 
that all sources with mass-loss rates at or below this value 
will have stopped expanding by this time. According to 
the expansion model, if the bubble remains radiative, at 
0.22 Myr the shell radius will be Rg = 0.63 pc. For a 
source with a ttIw = 2 x 10 “^ M0yr“^, the highest mass- 
loss source in W2, at the time of pressure confinement 
its size will be Rg = 0.88 pc. This is comparable to the 
radii that we find in the simulation, and Figure shows 
the total combined effective shell volume of the sources 
in W 2 is around Ipc at 0.2 Myr. Since runs W1 and W 2 
contain multiple sources, the net expansion is dominated 
by the strongest wind progenitor. As individual bubbles 
become pressure confined, the net expansion slows. The 
strongest sources have a mass-loss rate of 1.8 x 10 “^ and 
2.4 X lO“^M 0 yr“^, respectively, which correspond to a 
expansion time of 0.33Myr and 0.91Myr. The exact time 
to pressure equilibrium depends on the source location, 
i.e., variation of the ambient density and magnetic field 
in the star vicinity. 


3.2. Turbulent Dissipation 

We use the global velocity dispersion to estimate the 
impact of wind driving on the turbulent dissipation. Vi¬ 
sually, the slope of the velocity dispersion (for runs with¬ 
out external driving) in the top panel of Figure 0 indi¬ 
cates the rate of dissipation. We note that the 3iange 
in the global velocity dispersion for the run with neither 
external driving nor winds following the turbulence ini¬ 
tialization (W0_T2) is very similar to those with winds 
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TABLE 2 
Bulk Results^ 


Model 


Pm 

fv 

/p,i,CTl£) > 2.5 


fE 


W1_T1_L0 

1.12 

1.23 

0.11 

0.011 

0.083 

0.082 

0.39 

W1_T1 

1.13 

1.24 

0.11 

0.011 

0.082 

0.082 

0.40 

W0_T2 

1.00 

0.92 

0.00 

0.018 

0.0079 

0.037 

0.016 

W1_T2 

1.05 

1.25 

0.094 

0.018 

0.094 

0.037 

0.44 

W2_T2 

1.00 

1.00 

0.016 

0.018 

0.018 

0.037 

0.063 

W2_T3 

1.00 

1.00 

0.015 

0.014 

0.017 

0.029 

0.060 

W2_T4 

1.00 

0.95 

0.013 

0.020 

0.015 

0.054 

0.058 

W3_T2 

1.00 

1.00 

1.5 X 10-^ 

0.018 

0.018 

0.037 

0.037 

W4_T2 

1.00 

1.00 

2.1 X 10-4 

0.018 

0.018 

0.037 

0.037 


^ Model name, the ratio of the rms magnetic field strength at tf =0.1 Myr to = 0 Myr, the ratio of 
the rms 3D Mach number at to volume filling fraction of the wind ai tf (defined as the volume 
of cells for which T > lOOK), momentum fraction of cells with a ID rms velocity > 2.5km s“^ at the 
initial and final times, and the energy fraction in cells with a ID rms velocity > 2.5km s“^ at the initial 
and final times. 


but no external driving (W2_T4, W1_T2, and W1_T1). 
More quantitatively, we define the average turbulent dis¬ 
sipation simply as Eturb = different times. 

This simply reflects the global change in energy and is an 
indirect measure of the dissipation which occurs on small 
scales. Run W0_T2, which has no winds or turbulent 
energy injection, provides an estimate of the expected 
dissipation rate, Eturb — 6 x 10^^ erg s“^. The total ini¬ 
tial turbulent energy is Eturb = 1-47 x 10^7 erg, which 
gives a dissipation time of tdiss = Eturb/Eturb ^ 0-8 
Myr. If the dissipat ion time is tdiss = ( [McKee 


19891 [Mac Low 1999|), then 77 0.5, which is within 

the expected range. Here, the free-fall time is given by 


te = y(37r/(32Gp)) =1.47 Myr. 

The model W1 velocity dispersion is initially enhanced 
by the winds, but begins to show strong dissipation after 
0.05 Myr. Using the times between 0.1 Myr and 0.15 
Myr, we estimate a dissipation rate of ^turb = 1.68 x 10^^ 
erg s“^, which is twice as high as the non-wind case. If 
the dissipation rate scales with the rms velocitv, E oc 
( [Mac Lo^|l999| ), the dissipation rate in run W1 at 
O.lMyr should indeed be (3crvui/crvE2)^ = (2.5x 10^/2.Ox 
105)3 ^ 2 times larger than run W0_T2 at the same time. 

The energy injection rate from all winds is Eyj = 
^'Eimuj^iv‘^ Although the wind velocity varies as a 
function of mass, models W1 and W2 have maximum 
wind velocities of 200 kms“^, such that the wind speed 
is effectively 200 kms“^ for all stars. Models W1 and W2 
have an energy injection rate of 5.6 x 10^^ erg s“^ and 
5.2 X 10^^ erg s“^, respectively (see Table [B. Both val¬ 
ues exceed the global energy dissipation in me non-wind 
simulation, which initially might imply that the wind in¬ 
put should be sufficient to maintain the global turbulent 
energy. However, since the wind expansion slows as the 
shells become pressure confined, the high wind velocities 
and energy input do not directly translate into turbulent 
motion. 

In Run W2_T4 the rms velocity goes up by ^ 0.1 
percent initially, but turbulent dissipation quickly dom¬ 
inates. From 0.02-0.1 Myr, the turbulent dissipation is 
5.2 X 10^^ erg/s. This is about 15% smaller than the fidu¬ 
cial value, however more strongly magnetized turbulence 
has a smal ler dissipation rate (e.g. Stone et al.|1998 Mac 


Low 1999), and this value is consistent with the prior 
numerical studies of turbulent dissipation. Although the 


winds do inject energy locally, these results suggest that 
they are not sufficient to offset the turbulence dissipation 
and cannot maintain the initial cloud turbulence. 


3.3. Turbulent Statistics 

In order to assess the impact of the winds on the 
cloud properties and distribution, we compute several 
common turbulent statistics: the Fourier velocity, den¬ 
sity, and momentum power spectra and the density, ve¬ 
locity and momentum probability distribution functions 
(PDF). The choice of these statistics is partially moti¬ 
vated by prior work indicating their utility as di agnos¬ 
tics for various turbulent properties. For example, [Swift 
& Welch (2008) found that bipolar protostellar outflows 


can generate a feature in the po wer spectrum of the 


CO i ntegrated intensity. Similarly, Carroll et al. (2009 
2010) demonstrated that outflows niay impact th e veloc 


ity power spectrum slope. Federrath et al. (2010) showed 
that the mode of turbulent forcing, whether compressive 
or solenoidal (stirring) also impacts the gas density PDF, 
where compressive forcing produces a broader range of 
densities. Here, we examine whether the wind shells pro¬ 
duce additional compressive forcing. 


3.3.1. Probability Distribution Functions 

The density PDF of gas in driven turbulence simula- 
tions is generally very close to a log-normal distribution 


([Mac Low|1999[ |Padoan fc N ordlund|2004| [Kritsuk et al. 

2007| LUderrath et al.[[20in): 

1 


PiCjdt, = 


^/2 


-exp 


TTCF^ 


2cr2 


df, 


(15) 


where ^ = ln(p/p), a is the width, and normalization 
ensures that J p{^)d^ = 1. The top panel of Figure 
shows the density PDFs with and without winds. The 
non-wind turbulent density PDFs resemble a log-normal 
distribution as expected. However, once wind launch¬ 
ing commences gas is evacuated from a fraction of the 
domain, leaving low-density voids. This manifests as a 
second peak at densities of p ^ 10“^^ — 10“^^ g cm“^. 
However, gas redistribution appears to have little impact 
at high-densities, i.e., the mass that piles up in the shells 
does not modify the global density distribution on large 
scales. 

The velocity PDF, shown in the middle panel of Figure 

exhibits the most difference in the case with winds. It 
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Fig. 3. — Log density slices through four star locations in run W1_T2 at O.lMyr. The vectors show the direction and relative strength 
of the magnetic field. The star locations are marked with white crosses. 


TABLE 3 
Spectral Slopes^ 


Model 

t(Myr) 

ap 


CXy 

W2_T3 

0.0 

-1.39 

-0.58 

-1.7 

W2_T3 

0.1 

-1.38 

-0.55 

-2.5 

W1_T1 

0.0 

-1.38 

-0.50 

-1.7 

W2_T2 

0.0 

-1.33 

-0.56 

-1.7 

W2_T2 

0.1 

-1.30 

-0.50 

-2.6 

W1_T2 

0.1 

-1.43 

-0.45 

-2.4 

W1_T1 

0.1 

-1.57 

-0.39 

-2.4 

W2_T4 

0.0 

-1.39 

-0.39 

-1.7 

W2_T4 

0.1 

-1.41 

-0.39 

-2.7 


^ Model name, time, slope of the Fourier 
spectrum of logarithmic density, slope 
of the Fourier spectrum of momentum, 
and slope of the Fourier spectrum of ve¬ 
locity. The analysis employs a uniform 
512^ grid, which is interpolated from the 
AMR data. The fits are performed over 
the range k = 4.3 — 16.7. The 1 — cr un¬ 
certainty in the ap, and ay fits are 
approximately ±0.02, ±0.05, and ±0.1, 
respectively. We list the runs in the ta¬ 
ble in order of increasing magnetic field 
(see Figure^. 


develops a high-velocity tail that extends to several hun¬ 
dred knis“^. The peak at 200 knis“^ corresponds to the 
wind launching velocity maximum set in the simulations. 
As expected, the runs with the highest total mass-loss 


generate the largest amount of fast moving, low density 
gas. Otherwise, the distributions remain very similar to 
the initial turbulent PDFs. 

The third panel of Figure shows the momentum 
PDF. All momenta PDFs appear to be approximately 
log-normal, with similar widths. Some asymmetry in the 
density and velocity distribution appears as a result of 
the random driving. The PDF of the gas with and with¬ 
out winds is nearly indistinguishable. Essentially, the 
wind shells are shocks, which have very similar momenta 
to the initial turbulence. The shell expansion velocities 
are ^1 — 3 kms“^, which bracket the turbulent veloc¬ 
ity peak at 1.2 kms“^. This accounts for the minimal 
change in the bulk velocity and momentum distributions. 
Likewise, the PDFs for the strong magnetic field case are 
not statistically different than the weaker field case. 

If the winds are preferentially adding compressive 
modes, then we might expect the density and momenta 


PDFs to significantly broaden (e.g., Federrath et al. 2008 


2010). However, aside from the appearance of a seconc 


peak in the density PDF, this does not occur to any sig¬ 
nificant degree. Quantitatively, the standard deviation of 
the momentum distribution increases by 7% and 0.8%, 
respectively, for runs W1_T1 and W2_T2. The outputs 
with winds are otherwise similar to the pre-wind turbu¬ 
lence, even in the case with the highest mass-loss rates, 
in which the winds have a volume filling fraction > 10%. 
The similarity of the pre- and post-wind runs could also 
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0,00 0,05 0,10 0,15 0,20 

time (Myr) 


Fig. 4.— Velocity dispersion (top), rms magnetic field (middle) 
and effective radius ( defined as i^eff = (3V/47r)^/^ for the total 
volume enclosed by shells, bottom) as a function of time for the 
various models, where Brms = Si(B^ ■ + - and 

i^eff = [^*^^t>ioook/(^^/^)] legend also displays the 

total mass-loss rate for each run in MQyr“^. The vertical black 
dashed lines indicate the time that a source with mass-loss rate 
Mw = 10“®, 10“^MQyr“^, respectively, is predicted to become 
pressure confined for our fiducial parameters (B = 13.5/rG, Cg = 
0.19kms“^, (JiD = 1.14kms“^, and Vw = 200kms“^). 


be b ecause the winds do not create non-local effects (see 


3.3.2. Power Spectra 

We compute the Fourier power spectrum of t he density, 
mom entum and energy distribution using (e.g.|Federrath| 

P{q,k) = {q-q*^T^k'^)k^ ( 16 ) 

where q is the Fourier transform of a three-dimensional 
quantity q{l): 


q(k) = 


1 


(27rL)3/2 




(17) 


for I in the range [0, L] with the maximum spatial scale 
L corresponding to the smallest wavenumber k = 27rlL. 

Although the effective grid resolution including AMR 
is 1024^, the basegrid resolution (minimum resolution 
without AMR) is 256^. This means the initial inertial 
range is /c ^ 5 — 20, since there is little AMR refinement 
at early times. The inertial range is limited at small 
k because we inject turbulence over the wave numbers 
k = 1 — 2. Once wind driving begins, the wind shells 
and strong shocks will be refined to higher resolution, 
increasing the effective inertial range. We perform the 
analysis using flattened 512^ grids for all the variables. 
A comparison of 256^ and 512^ resolution is presented in 
Appendix We summarize the best-fit slopes for each 
of the Fourier spectra in Table 

The top panel of Figure shows the Fourier spec¬ 
trum of logarithmic density. The spectral slopes are 
comparable to S{k) oc with the winds causing a 

slight steepening in the slope in the cases with strong 
winds (Wl). The fitted slope s are all consistent with 
the solenoidal driving case in Federrath et al. (2010), 
for which ap oc k~^‘^. Compressively driven turbu¬ 
lence exhibits a significantly steeper density spectrum: 
ap oc k~‘^-^. Since the W2 snapshots with winds are not 
significantly statistically different over the inertial range, 
this suggests the turbulence is still largely solenoidal. 
The slope of Run Wl does steepen slightly, which may 
indicate moderately enhanced compression. The winds 
impact the shape of the spectrum over the dissipation 
range (i.e., k = 40 — 100), which becomes slightly flatter 
in the cases with winds. This does not appear to be a 
product of the increased resolution in the wind runs since 
the same difference appears for power spectra computed 
at a lower grid resolution (see Appendix IS and for the 
run W1_T1_L0. 

The Fourier spectra of the momentum (g = pv) are 
shown in the middle panel of Figure These spectra are 
likewise very similar over the inertial range for the wind 
and non-wind cases. The only difference attributable to 
winds manifests over the dissipation range at larger k val¬ 
ues. The outputs with winds are again slightly system¬ 
atically higher (more momentum on small scales) than 
the same model without winds. 

The velocity power spectra (bottom panel of Figure 
is the only statistic to show a strong departure from the 
non-wind case. A peak occurs at low k values, /c ^ 4 —10, 
which is indicative of the bubble size, with the strong 
wind case peaking at the smallest k (/ ^ A/4). This 
creates a pronounced steeping of the slope, ay —2.4, 
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Fig. 5. — Density (left), velocity (middle) and momentum (right) PDFs for different models at times prior to wind launching (0 Myr, 
lines) and at 0.1 Myr (lines and symbols). A lognormal function in grey is shown in the left plot: p(lnp) = 9.7 x lO^exp [—1.3(lnp — 21.2)^]. 


for a power-law fit over the inertial range. As in the other 
cases, the winds impact the dissipation scale by adding 
a little extra power on small scales. 

The power-law slopes are nearly independent of the 
average magnetic field strength. Only the Fourier mo¬ 
mentum spectral slope exhibits a statistically significant 
change: the slope flattens with increasing magnetization. 
However, the wind feature in the velocity spectrum no¬ 
ticeably shifts to smaller k with decreasing field strength. 
This is indicative of the role of magnetic pressure confine¬ 
ment in setting the bubble size. Consequently, stronger 
field strengths correlate to slightly steeper slopes. 

All the velocity slopes are flatter tha n the solenoidal 
and co mpressive driving slopes found by [Federrath et alT 
(2010), ay -1.8 to -2.0 for non-magnetized turbulence. 


In general, the veloci ty power spectral s lope is ex pected 
to fall between — | (Kolmogorov 1941) and —2 (Burg- 
ers||1948|). The difference may be partially due to lower 
resolution, which is discussed in Appendix [B] 

3.4. Synthetic CO Distributions 
3.4.1. CO Modeling Method 

We post-process the simulations with the radiative 
transfer code radmcSeQ to obtain the emission in the 
CO(1-0) molecular line. We compute the emission us- 
i ng the Large Velo city Gradient (LVG) approximation 
(jShetty et al | 2011 |), which calculates the rotational level 
populations by solving the equations for local radiative 
statistical equilibrium. We perform the radiative transfer 
using the densities and velocities on a 256^ grid, where 
we define = p/( 2 . 8 m^) and adopt a C O abundance 
of [^^C 0 /H 2 ] =10“^ ( [Frerking et al. 1982). Gas above 


800 K or with < 10 cm“^, i.e. effectively all gas 
in the rarefied bubble regions, is assumed to have a CO 
abundance of 0 . Gas with densities uhs > 2 x 10 ^ cm~^ 
exper iences freeze-out onto dust grains (Tafalla et al. 
20041) , which significantly reduces the gas phase CO. We 
set tne abundances to zero in these high density regions. 
We include turbulent line broadening on scales at and 


^ http://www.ita.uni-heidelberg.de/ dullemond/software/radmc- 
3d/ 


below the grid resolution by adding a constant micro¬ 
turbulence of 0.25km s“^, which is dete rmined by the 
observed linewidth-size relation (see 12.2). The spectral 
resolution of the output cube is At; = 0.156 kms“^. 

To model an observation of Perseus, we assume the 
simulated clouds are located at a distance of 250 pc and 
convolve the emission with an 46” beam, which is the res- 
oluti on of the COMPLET E ^^CO (1-0) survey of Perseus 
(e.g., |Ridge et al.||2006b|). 

Since the morphology of the observed bubbles suggests 
that Perseus is flatter along the line of sight than its 
projected length (e.g.. All conclude dos ^ 0.5pc), we 
also calculate a spectral cube for a thin cloud, which is 
8i Lz = 0.8 pc wide subset of the whole domain (Run 
W1_T1_S1). We choose the subdomain such that some 
of the bubbles span the entire cloud depth and some 
sight lines contain no molecular emission. Indeed, some 
of the larger bubbles observed in Perseus “break out” 
of the cloud and do not enclose molecular gas. Similar 
cloud-s panning bubbles have also been identified in other 
clouds jBeaumont fc Williams 2010 Debarveng et al. 


2010 INakamura et al.|| 2012 p . 


3.4.2. Synthetie CO Intensity Maps 

Figure shows maps of the integrated brightness tem¬ 
perature. Some of the wind shells appear clearly in the 
CO emission, although not all resemble rings due to su¬ 
perposition or interaction between the winds. The CO 
line is strongly excited in the shells because the gas is 
both denser and warmer than the mean cloud density 
(nH 2 — 400 cm“^) and temperature (T = 10 K). At 0.1 
Myr the shell sizes are typically less than 1.5 pc with a 
strong dependence on the source mass-loss rate and some 
dependence on the ambient density near the source. The 
projected radii of the CO shell emission in Perseus are 
^ 0 . 2 - 2.8 pc, which is comparable to the sizes of the shells 
modeled here. Given an expansion rate of 0.1 pc per 0.1 
Myr for the strongest source in model W2 (see Figure^, 
the radius of the largest shell might exceed 2 pc following 
another Myr of evolution. 

Although the shells do not resemble circular rings, 
they are distinctly rounded structures. The shells ap- 
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Fig. 6. — Log density (top), momentum (middle) and velocity 
(bottom) Fourier spectra computed using a 512^ grid resolution. 
The grey line denotes the approximate inertial range and indicates 
the typical best fit slope. The density and energy spectra are com¬ 
pensated by kP ‘. The input quantities are normalized to the average 
density, rms momentum, and rms velocity, respectively. 


pear more clearly in the case where the bubbles span 
the width of the cloud. Like observed shells, the simu¬ 
lated structures have variable brightness, i.e., they are 
clumpy. Here, this is caused by density variance in the 
turbulent gas, such that different parts of the shell in¬ 
teract with denser ambient material. If the ambient gas 
is lower density then the hot gas expands more freely, 
leading to additional asymmetry. These characteristics 
agree with the observations of Perseus, for which not 


all shells are complete circles or are easily identified in 
the integrated emission map; some must be identified on 
the basis of their coherent velocity structure in the spec¬ 
tral cube. The simulated sources with weaker winds, 
< lO“^M 0 yr“^, also produce much weaker shell fea¬ 
tures, which are not necessarily visable in the CO emis¬ 
sion. 

Wind model W 2 , which has sources with similar mass- 
loss rates as those estimated in Perseus, also exhibits 
shell structure with similar CO intensity. As reported in 
All, the identified Perseus shells have average CO inte¬ 
grated intensities of ^ 0.6 — 2.4 K kms“^ (see their Fig¬ 
ure 4). The shells in the simulated W2 maps, which have 
been convolved to similar spatial and velocity resolution 
as the observations, have CO integrated intensities in the 
range of 0.7 to 2.1 Kkms“^, similar to that observed in 
Perseus. Without winds, the CO integrated intensities 
do not exceed 0.7 K kms“^. In the strong wind case, 
Wl, shell emission ranges from 0.7 —7.0 K kms“^, which 
exceeds the shell emission in Perseus. This supports the 
estimate by All that sources in Perseus have mass-loss 
rates from lO“^M 0 yr“^ to a few times lO“^M 0 yr“^, 
which is consistent with run W2. For wind velocities 
of 200km s“^, model W2 has an input momentum of 
P = Mvw = 9 X lO“^M 0 kms“^yr“^. Although we can’t 
directly confirm that all obser ved momentum is gener¬ 
ated by line-driven winds (see ^4.3), we have high confi¬ 
dence in the magnitude of the momentum input needed 
to produce the observed features. The stronger sources 
in model Wl produce features that are too bright, so it 
is unlikely that any of the Perseus sources have mass-loss 
rates exceeding 5 x lO“^M 0 yr“^, in good agreement with 
the observed estimates. 


3.4.3. Hot Gas Distribution 

As discussed by All, some of the identified shells are 
ambiguous. Either their signature in the velocity chan¬ 
nels is less clear or there is no obvious co-spatial wind 
progenitor. Turbulence, which naturally produces spa¬ 
tial fluctuaions over large scales, provides one alterna¬ 
tive explanation for these features. Statistical fluctua¬ 
tions in the density and velocity distribution may pro¬ 
duce voids, which could masquerade as wind features. 
To test this possibility, we examine the distribution of 
hot gas. Figure shows an integrated mask selecting gas 
with temperatures greater than 800 K. Since gas is only 
heated through wind interactions, gas temperature serves 
as more direct tracer of shells than gas density. Compari¬ 
son of Figures and show that most sources are indeed 
associated with regions of hot gas. However, in many 
cases the hot gas is able to leak out of the immediate 
vicinity of the source and create significantly extended 
and asymmetric bubbles. The hot gas quickly expands 
in lower density regions, so the bubble shapes are in some 
sense indicative of the pre-wind turbulent structure. The 
asymmetric expansion of the bubbles leads to the appear¬ 
ance of distinct shells without any central source (e.g., 
top-right feature in the images showing run W1_T1). Be¬ 
cause the domain is periodic some hot gas also appears 
on the opposite side of the domain, which highlights the 
large extent of some of the bubbles. 

Comparing the hot gas distribution in Firare with 
the apparent voids (dark regions) in Figure indicates 
that many but not all apparent voids are associated with 
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hot wind gas. For example, VI and V2 indicate two 
evacuated regions that do not have corresponding hot 
gas in Figure This suggests that some small fraction 
of the shells may be caused by turbulent stochasticity 
rather than winds. Unfortunately, the temperature of 
low-density material is difficult to probe observationally 
when there is little or no CO emission present. Further¬ 
more, even in cases where atomic emission is detected, 
the gas may be warm because it has been dissociated by 
the background stellar FUV field rather than by heated 
by local sources. 


3.4.4. CO Distribution 

Figure shows the rms velocity distribution of the 
smoothed CO emission for snapshots with and without 
winds. The rms velocity, 'r’rms: is defined as the second 
moment of the intensity distribution along a given sight- 
line: 


"^rms — 


-Vcent)‘^/ 


0.5 


(18) 


where 'r’cent is the centroid velocity or first velocity mo¬ 
ment: 

'^cent ~ (19) 

As illustrated by Figure winds create a distinctive tail 
in the velocity distribution. Although most of the tail 
is only a few kms“^, it extends the distribution signif¬ 
icantly above the median cloud velocity. This range of 
velocities, 3 — 5 kms“^, is also consistent with the mea¬ 
sured expansion velocities of observed shells (e.g. All). 
We also note that the shape of the velocity feature is 
similar t o that of the observed velocity distribution in 
Perseus (Beaumont et al. 2013), where it was conjec¬ 
tured to arise from stellar feedback. Table gives the 
mean, median, and standard deviation of the v^ms distri¬ 
bution. The cases with winds have systematically higher 
mean and median Vj-ms values and have standard devia¬ 
tions as much as three times higher. We note the highest 
velocity material generated by the winds occurs in hot 
(ionized or atomic) gas and thus is assumed to have zero 
CO abundance by construction. 

Assuming that the integrated intensity is a proxy for 
mass, we can estimate the fraction of momentum and 
energy in high-velocity gas. Prior to wind launching less 
than 0.1% of the momentum or energy is in emitting gas 
with a line-of-sight velocity Vj-ms > 1km s“^. However, 
Table shows that winds introduce higher velocity ma¬ 
terial and a significant fraction of both the momentum 
and energy falls above this threshold. The exact amount 
of high-velocity gas is very sensitive to the total mass-loss 
rate, with up to 50% of the momentum being in winds for 
the highest mass-loss rates. For a mass-loss rate similar 
to that estimated for Perseus, 10% of the momentum is 
above this threshold and 20-30% of the energy. The cloud 
slice (W1_T1_S1) and full cloud have a similar fraction 
of energy and momentum above the threshold. 


3.4.5. CO Turbulent Statistics 

Occasionally, analysis of the integrated CO emission 
or column density can provide insight into t he under- 
lying turbulence or star formation activity (Swift & 


Welch 2008 Carroll et al. 2010 Burkhart et al. 2010). 
NGClo33, wnich exhibits copious outflow activity shows 
a spectral break around /c 20 in CO(1-0). No break 
is apparent in the ^^CO emission, although this may be 
attributed t o the higher optical d epth of ^^CO. 

Following Swift & Welch (2008) we calculate the power 
spectra of the integrated intensity for ^^CO and ^^CO. 
This entails performing a Fourier transform of the 2D 
intensity map, squaring the result, and then averaging 
over annular rings in k space. We subdivide the in¬ 
tegrated emission into red-shifted and blue-shifted gas, 
which spans the range of shell velocities. For the simu¬ 
lated clouds, there is no net global velocity, so the red- 
and blue-shifted ranges are simply defined symmetrically 
about 0 kms“^. 

Figure shows the power spectra of this emission for 
various wind models. We find the red- and blue-shifted 
emission have similar slopes and none of the models with 
winds show any clear feature in the spectrum that indi¬ 
cates a preferred energy injection scale. The ^^CO line 
wings of the model with only turbulence (t = OMyr) have 
slopes of —2.7 ± 0.1 and —2.8 ± 0.1, while W2_T2 and 
W2_T3 have slopes ranging from —2.5 to —2.7. These 
slopes are slightly shallower, but consistent within la er¬ 
ror. The strongest wind case, Wl, has ^^CO line wing 
slopes of —2.3±0.2 and —2.54±0.05, which is marginally 
flatter than the weaker wind runs and could indicate a 
trend with increasing wind activity. 

The ^^CO line wings each have a slope of —2.35 ±0.02. 
This result is consistent with the expectation that opti¬ 
cally thick lines saturate at a slope of —3, w hile optically 


thiner lines, like have flatter slopes (Lazarian & 

Pogosyan|2004| iBurkhart et al.|2013|). Additional analy¬ 
sis exploring difierent velocity ranges, the full integrated 


intensity cube, and the gas column density likewise re¬ 
veal no clear break in the intensity spectrum. The corre¬ 
sponding slopes for different velocity ranges tend to vary 
by an additional ±0.1 with higher velocity channels with 
less average emission producing flatter slopes. 

While the lack of a feedback signature may be sur¬ 
prising, breaks in observed spectra are rare and seem to 
appear only in regions with very vigorous feedback, such 
as NGC1333. No observational data to date has revealed 
any such feature attributable to feedback from winds. 


4. DISCUSSION 
4.1. Wind Signatures 

Stellar winds create a significant amount of high- 
velocity gas that leaves a clear signature in the density 
and velocity PDF, however most of this gas is quite rar¬ 
efied and is confined to the central regions of the shells. 
This material, given its density and temperature, is not 
in molecules and would not show up in CO emission (ex¬ 
cept as a void). It might be possible to detect these 
regions in radio or X-rays, however. The cold molecu¬ 
lar material entrained and accelerated by the expanding 
shells does produce a clear signature, one that is consis¬ 
tent with the characteristics of the observed CO emission 
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Fig. 7. — Integrated ^^CO(l-O) brightness temperature for Run W2_T2 (top left), Run W2_T3 (top right), Run W1_T1 (bottom left) 
and a O.Spc slice of Run W1_T1 (W_1T1_S1, bottom right) at 0.1 Myr. The projected locations of the stars are marked with asterisks. VI 
and V2 mark two voids discussed in the text. 


TABLE 4 

Synthetic Emission Properties 


ModeE 

t(Myr) 

t’(kms ^) 

i;med(kms 0 

Gv (km s ^) 

^Erms|>l/^tot 

^Erms|>l/^tot 

W1_T1 

0 

1.41 

1.40 

0.29 

0.0 

0.0 

Wl.Tl 

0.1 

2.04 

1.77 

0.86 

0.51 

0.66 

W1_T1_S1 

0 

0.73 

0.69 

0.26 

0.0 

0.0 

W1_T1_S1 

0.1 

1.02 

0.82 

0.73 

0.29 

0.50 

W1_T2 

0.1 

2.06 

1.67 

0.97 

0.63 

0.77 

W2_T2 

0 

1.43 

1.41 

0.29 

0.0 

0.0 

W2_T2 

0.1 

1.49 

1.43 

0.44 

0.10 

0.21 

W2_T3 

0 

1.33 

1.33 

0.25 

0.0 

0.0 

W2_T3 

0.1 

1.45 

1.35 

0.53 

0.14 

0.30 

W2_T4 

0 

1.43 

1.35 

0.30 

0.0 

0.0 

W2_T4 

0.1 

1.46 

1.39 

0.48 

0.11 

0.25 


^ Model name, mean rms spectral velocity, median rms spectral velocity, standard deviation of the 
rms spectral velocity, fraction of momentum in pixels with los rms velocities |nrms| > 1km s“^, and 
fraction of energy in pixels with los rms velocities |i;rms| > 1kms“^. We assume the CO intensity 
is a proxy for the underlying gas mass i.e., we assume the emission is optically thin, which is the 
assumption that underlies the CO X-factor. 

in Perseus. 

There are other clear differences between the bulk 
properties of the CO emission and the underlying den¬ 
sities and velocities. For example, the total amount of 
momentum and energy added by the winds is only a few 
percent (except in the highest mass-loss case). Mean¬ 
while, the total energy and momentum inferred to be in 
the high velocity gas based on the CO emission is 20- 
30%. As much as 50-60% of the momentum and energy 
appears to be in the shells in the highest mass-loss run. 


This value is in fact comparable to the total fraction of 
cloud energy inferred to be in these shells for Perseus 
(50%). 

The discrepancy between the information obtained 
from the CO emission and 3D kinematic information. 
Tables]^ andrespectively, underscores that CO is a bi¬ 
ased probe ofthe total underlying cloud momentum and 
energy. It preferentially selects optically thin gas with 
densities of nHs — 500 — 5 x 10^ cm“^, and the aver¬ 
age wind shell densities fall within this range. The shells 
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Fig. 8 .— Integrated temperature mask including all gas with T > 800 K for Run W2_T2 (top left), Run W2_T3 (top right), Run W1_T1 
(bottom left) and a 0.8pc slice of Run W1_T1 (W_1T1_S1, bottom right) at 0.1 Myr. The colorscale indicates the fraction of the cloud 
along each sightline with temperature above 800 K. The projected locations of the stars are marked with asterisks. The image is saturated 
where the shells have broken-out of the cloud. 


are also slightly warmer than the exterior gas, which en¬ 
hances the CO excitation. We conclude that while CO 
emission is a useful diagnostic of cloud structure and may 
be used to identify the location and properties of shells, 
it likely over-estimates their impact compared to the 3D 
gas densities and velocities. 

The similarity between the simulated and observed CO 
emission leads to several practical conclusions. First, the 
mass-loss rates inferred from the observations are likely 
accurate; without energy injection of this magnitude it 
would be impossible to create the observed shell sizes and 
energies. Second, these mass-loss rates also reproduce 
the energy and momentum inferred from observed CO 
emission, which conhrms that shell progenitors do play a 
signihcant role in shaping the gas morphology and ener¬ 
getics. Finally, since winds produce a clear signature in 
the molecular emission, evidence of this feature in obser¬ 
vations is indicative of the presence of stellar winds. 


4.2. Turbulent Driving 

A fundamental question motivating this study is “Do 
stellar winds drive turbulence?” It is clear that signih¬ 
cant gas motions do result as a product of the winds, 
but it is less clear whether these motions are purely lo¬ 
cal and, consequently, do not contribute much to the 
global turbulence. For compressible turbulent huids, en¬ 


ergy cascades fr om larger to smaller scales (e.g. Sridhar & 
Goldreich||1994 ), where the largest energy injection scale 


is the peak of the cascade. Little or no energy propa¬ 


gates from smaller to larger scales (a “reverse” cascade). 
In the case of protostehar outhows, shocks created by 
the bipolar hows have a strong local impact, changing 
the fraction of solenoidal motions of the turbulence and 
creating new motions on small scales below the size of the 
outhow ene rgy injection scale (e .g., the linear extent of 
the outhow Ohner & Arce|2014|). Simulations show that 
within parsec sized clusters outflows can contribute suffi¬ 


cient turbulence to halt global collapse and regulate the 


star formation rate ( 

Nakamura & Li||2007 

1 Carroll et al. 

20091 Wang et al.|20] 

Lo Hansen et al.|20I2[ 

). T'he dynam- 


free fall time and it appears that outhows can replenish 
turbulence on scales that impact incipient star forma¬ 
tion. Likewise ionization from high-mass stars has a sig¬ 
nihcant impact on cloud structure and appears to be able 


to inject su fficient energy to drive turbulence (Boneberg 
et al.|[2M5 ). In contrast, the winds are largely spherical 


and sweep up gas as they advance. However, two ehects 
limit the efficacy of energy injection by winds. First, 
while the energy injection scale, i.e., the ehective radius 
of the bubble, increases with time, the bubble expan¬ 
sion stalls for mass-loss rates of < yr“^. This 

ehect is enhanced by stronger magnetic helds. Second, 
high-velocity motions that are created are in the low- 
density material behind the shock, which remains inside 
the shells. This residual gas contains very little mass and 
hence does not contribute much to the global turbulent 
cascade. 



















15 




(j 

CL 



Fig. 9.— Histograms of the intensity-weighted rms velocity, i;rms, 
for all sight lines through the cube. The spectral cubes have been 
convolved with a beam of 46” at d=250pc. The symbols indicate 
the snapshots with winds. 


Whether winds can drive turbulence non-locally de¬ 
pends upon whether there is any impact on gas motions 
outside the shell radius. Alfven waves may provide one 


avenue to excite non-local motions (e.g., Heitsch & Burk- 


ert 


2002). The expanding shells displace the magnetic 


held, which excites Alfven (transverse) waves that prop 
agate along the field lines with a velocity va = B j 
where B is the magnetic field strength and p is the 
gas density. For the average simulation field strength 
and gas density, Alfven waves propagate with a veloc¬ 
ity of 0.8 kms“^, which is less than the shell expansion 
velocity of 1-5 kms“^. Fast magnetosonic waves may 


also transport ener gy away from the shells (Gendelev 


& Krumholz 2012 ). Fast waves propagate perpendic- 


ular to the field with a speed Vf = \/v\ + , which is 

faster than the Alfven speed but not significantly so since 
the sound speed in the gas just outside of the shells is 
Cs — 0.2 kms“^. Due to the turbulent nature of the 
gas and tangled field morphology, it is not possible to 
di rectly see the propagation of magnetosonic waves as 


Gendelev & Krumholz ( 2012 [ ). Figure 11 shows the 


m 

distribution of Alven velocities with and without wind 
driving. The shocked shell gas (p > 10“^^ g cm“^) 
has Alfven velocities smaller than the shell expansion 
speed. Any excited magnetic w aves may also rapidly 
damp ( Banerjee fc Pudritz|2007 ). Consequently, we find 
the runs with winds only exhibit a minor enhancement 
in field strength and Alfven speed at high densities. Lit¬ 
tle change occurs at intermediate and lower densities. 
The evacuated regions have densities < 10“^^ g cm“^ 
and a relatively small magnetic field. In summary, this 
suggests that shells do not contribute to the turbulent 
cascade except very locally. 


4.3. Scenarios for the Shell Origins 

Our simulations confirm the inferred mass-loss rates 
and wind momenta required to produce observed shells 
within low-mass star-forming molecular clouds. How¬ 
ever, these results reinforce the severity of the weak wind 


Fig. 10.— Power spectrum of the integrated CO intensity for 
gas with velocities —4.5 < v < 1.0km s“^ (blue-shifted) and 1.0 < 
V < 4.5kms“^ (red-shifted). The blue-shifted spectrums are offset 
by a factor of 50 for clarity. Power spectra are calculated for both 
^^CO(l-O) and ^^CO(l-O) (green-dashed line) for W2_T2 at 0.1 
Myr. 


problem and pose a significant challenge to our under¬ 
standing of stellar feedback. Here, we compare and eval¬ 
uate explanations for these features. 


4.3.1. Line-Driven Winds 

The primary explanation for these features, line-driven 
winds from main-sequence stars, is marked by a primary 
physical limitation. Photon driven winds are limited by 
the amount of momentum in the radiation field, which 
is set by the stellar luminosity, L. If all the photons 
radiated by the star are absorbed in the wind then the 
momentum injection rate is M^Vw = Lje. This assumes 
that only the first scattered photon contributes to the ra¬ 
diation pressure (“the single scattering upper limit”). If 
multiple photon scatterings contribute then the momen¬ 
tum injection may be boosted by a factor of the optical 
depth of the wind, where Ty^ is simply assumed to be 
constant with frequency. Thus, if the observed shells are 
created predominately by line-driven stellar mass-lo ss, 
then the mass-loss rate is ( Lamers fc Cassinelli||1999 ): 




TwL 


VwC 
= 2 


X 


X lO-^Moyr-i (^) 


L 


lO^L^ 


V 200 km S' 


( 20 ) 


where 10^ Lq is approximately the luminosity of a BO.5 
star and the lines are assumed to be optically thick. This 
estimation is sufficient to account for the inferred mass- 
loss rate for the strongest source {ttIw = 2 x lO“^M 0 yr“^) 
in both the W 2 simulation and in Perseus. 

In practice, main-sequence winds are driven in nar¬ 
row frequency bands, over which various atomic lines are 
optically thick. Higher luminosity stars have more sat¬ 
urated lines and hence can drive stronger winds, while 
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dimmer stars have fewer or no optically thick lines for 
which momentum can be deposited. This is the essence of 
the weak-wind problem. Thus, the above estimate, which 
requires high-optical depth and multiple-scattering, is 
optimistic. These complications suggest that another 
factor or additional mechanism likely enhances momen¬ 
tum deposition in this classic line-driven wind picture. 

4.3.2. Wind Variability 

While the average wind velocities and stellar mass-loss 
may be orders of magnitude lower than we predict here, 
it is possible that stars, especially very young ones, expe¬ 
rience periodic wind enhancements due to elevated mag¬ 
netic stellar activity or accretion. In fact, winds cannot 
be continuously three orders of magnitude higher than 
expected throughout the star lifetime since this would 
result in a significant change to the total stellar mass, 
e.g. of order ^ Mq. However, given a mass-loss rate of 
> IO^^Mq yr“^, the simulations demonstrate that the 
shells are created relatively rapidly. This suggests that 
the observed features could plausibly be caused by inter¬ 
mittent high-activity or accretion over 0.1 Myr. The 
cumulative affect of periodic wind enhancement could ac¬ 
count for the shell sizes where a prolonged epoch of low 
mass-loss alone would be insufficient. 

Magnetic activity or accretion enhancements may be 
caused by interactions with a close stellar companion or 
interaction with some nearby residu al gas or both (e.g., 
Basri et ar]|1997 ; Jensen et al.||2QQ7 ). Since B-type stars 


do not have a convection zone in their outer atmosphere 
they are generally not expected to be magnetically active. 
However, insta nces of magnetical l y active B-stars hav e 


been observed (Gagne et al.||I997 

Pillitteri et al.||20I4) 

A sub-type of B stars, /t-Cep varia 

Dies, do experience A 


ray pu lsations, which may a l so lea d to a wind enhance¬ 


ment (Oskinova et al. 2014 
which has produced a shel 


2015|). B-star p Ophiuchi, 
visible in the surrounding 


dust emission, is a bi nary and exhibits magnetic activity 
( Pillitteri et al.||2Q14|) . 

I'he multiplicity of the sources in Perseus is currently 
unknown. Stati stically, AB stars have a multiplicity frac¬ 
tion of 50-60% (Duchene & Kraus||2013), which suggests 
that at least hall of the wind progenitors may be in a bi¬ 
nary system. However, the separation distribution peaks 
around 50 AU, so the fraction of wind progenitors with 
close or eccentric companions (< lOAU) is considerably 
smaller. 

The gas distribution within the shells is also difficult to 
probe directly, since the gas is very low density and likely 
predominantly atomic. HI and H 2 data does not exist 
at the resolution required to probe the intra-shell gas. 
Contamination from the foreground and background ISM 
also make distinguishing gas associated with the shell 
features difficult. However, further observations of the 
source environments and properties could be revealing. 
All these possibilities warrant follow-up studies of the 
wind progenitors in Perseus. 

4.3.3. Protostellar Outflows 

A significant amount of energy and momentum is de¬ 
posited by magnetically collimated bipolar o utflows dur¬ 
ing the f irst ^ 0.5 Myr of the protostar stage. |Matzner~fc 
McKee (2000) estimate that the typical outflow momen¬ 


tum per unit stellar mass is Vc = 40km s“^. This is an ap¬ 
propriate value for a region like Perseus which, given the 
number of observed protostars and a standar d IMF, is 


likely not forming stars with masses above 3Mq (Matzner 
McKee fc Offner|20I0 ). The momentum deposition 


2007 


rate 


of the outflows can then be written 
dt 


m^Vc 


{tp) ’ 


( 21 ) 


where is the stellar mass and (tp) is the typical proto¬ 
stellar lifetime, over which bipolar outflows occur. The 
average rate of mass-loss due to outflows, rho, can be 
written 


m^Vc 

tpVo 


= 2.4 X 10“ 


—6 f 


( ) 

\40kms“^ / 


3 M 0 y \40kms“^7 \0.5Myr 


( __r' 

VlOOkms”^/ 


( 22 ) 

-1 

(23) 


where Vq is the outflow launching velocity. Alternatively, 
given some instantaneous accretion rate for a particular 
source, the X-wind and disk-wind outflow models predict 
that the outflow rate will be a fixed fraction of th e ac¬ 
cretion rate, fl, = 0.1 and fl, = 0.3 3, respectively (Shu 
et al.|I988t |Pelletier & Pudritz|I992|). Numerical models 
for disk wind launching hnd that the ratio of mass-loss to 
accretion rate can approach unity; however, there is an 
inverse correlation between mass loading and jet veloc¬ 
ity, so this may n ot correspond to an increa se in the total 
wind momentum ( Sheikhnezami et al.| 20 I 2 ). A protostar 
accreting at rh = I 0 “^MQyr“^ would thus be expected 
to have an outflow rate of rho — I — 3 x IO“^M 0 yr“^. 
The average accretion rate declines with time, so high- 
accretion is generally associated with the earliest stages 
of protostar formation, whe n the protostar is st ill em¬ 
bedded within a dense core (Dunham et al. |2QI4 ). 

While these mass-loss rates are comparable to those 
inferred for observed shells, the highly collimated nature 
of protostellar outflows requires that the output energy 
and momentum decline s teeply away from the outflow 
axis. iMatzner & McKee (2000) proposed that the force 


of the outflow deposition varies with angle according to 


p{d)c^ 


1 


ln(2/0o)[l + O 0 ~ cos 2 ( 6 ')] ’ 


(24) 


where 0 is measured from the outflow axis and Oq is the 
“flattening angle” of the force distribution, over which 
the bulk of the momentum is injected. Throug h com- 


parisons with observations, Matzner & McKee (2000) 


conclude that Oq < 0.05 and suggest a fiducial value of 
Oq = O.OI. These values indicate that the momentum ef¬ 
ficiency declines by a factor of 13-300 only 10° from the 
outflow axis. 

If we suppose that the shell progenitors have a bipo¬ 
lar outflow pointed largely along the line-of-sight, we can 
evaluate whether the transverse component of the out¬ 
flow would be sufficient to create the observed shells. 
For an outflow with an axis oriented 20° from the line- 
of-sight, any motions of the gas in the plane of the sky 
would be caused by outflow material 7(3° from the out¬ 
flow axis. At this angle assuming a conservative amount 
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of collimation, Oq = 0.05, the outflow force would be di¬ 
minished by a factor of fe = P(0°)/P(70°) ^ 350. For a 
SMq protostar, this corresponds to an effective momen¬ 
tum deposition rate of Pq = ferhoVo — 1/350 x (2.4 x 
lO“^M 0 yr“^)X 100 kms“^ ^ 5x lO“^M 0 yr“^kms“^ (see 
Equation [^. This amount of momentum deposition 
is significantly less than the momentum inject ion by 
a 3 M 0 star according to the Vink et al. (2001) mod¬ 
els: Pw = rh^Vw = lO“^M 0 yr“^ 6 OOkms“^ = 6 x 
lO“^M 0 yr“^kms“^. Runs W3_T2 and W4_T2, which 
modeled isolated stars with these properties, demon¬ 
strated that this momentum is too low to create struc¬ 
tures of the observed size and energy. 

In order for off-axis magnetically collimated outflows 
to create shell features, protostellar outflow rates would 
have to be tens to hundreds of times higher than those 
observed, i.e., ttIo ^ lO“^M 0 yr“^. Outflows of this mag¬ 
nitude are wildl y discrepant with inferred protostellar 
accre tion rates (Offner & McKee 2011 Dunham et al. 
2014), and there are no protostellar outflows iden tified 
within Perseus t hat have this amount of mass flux (Arce 
et al. 2010 2011). We conclude that protostellar outflows 
are very unlikely to produce the observed shells. 

4.3.4. Stellar Groups 

Small groups of stars may be able to provide sufficient 
integrated momentum to create a shell where a single 
star is insufficient. For example, Nakamura et al.|( 2012 |) 
conjecture that the shells near the LI641-IN and V 380 
Ori clusters were formed by the combined winds (or pro¬ 
tostellar outflow activity) of the cluster. Similarly, the 
shell surrounding rho Ophiuchi contains several addi¬ 
tional stars that may contribute to the total momentum 
(I. Pillitteri priv. comm.). However, this still presents an 
incomplete solution. The most massive star in L1641-N 
is a B4 V Herbig Ae/Be star, which has an extrapolated 
mass-loss rate < 1 O“^M 0 yr“^. The cluster contains 80 
young stellar objects, but these cumulatively are unlikely 
to account for the needed momentum from stellar winds, 


even extrapolating the Vink et al. (2001) prescription to 
lower masses. 

Although stellar groups may account for the properties 
of other observed shells, those in Perseus do not appear 
to contain star clusters. Shell progenitors are at most 
binary or triple systems. It is possible future surveys 
of stellar kinematics may reveal other associated stars, 
which were previously overlooked. 


4.3.5. Voids 

If voids are created by chance fluctuations in the turbu¬ 
lence then it is possible that low-density, coherent regions 
occasionally appear within clouds which are not due to 
feedback. If the gas is sufficiently low density then CO 
may not be present, i.e., the cloud is porous and the 
molecules are dissociated, or the emission may not be de- 
tectible. Thus, voids could masquerade as wind features. 
Our synthetic CO maps do show a couple of regions with¬ 
out CO emission that are not created by stellar winds. 
However, enhanced CO on the boundary of these voids 
is hard to produce without a central source. The sim¬ 
ulations don’t provide evidence that chance alignments 
are likely to create similarly bright morphological fea¬ 
tures. Consequently, we conclude that voids created by 



Fig. 11. — Average Alfven velocity (top) and average rms mag¬ 
netic field strength (bottom) as a function of density. The averages 
are computed as a density weighted average over each bin. The er¬ 
ror bars indicate Tier standard deviation in the binned values. For 
clarity, the dispersion is only shown for two representative cases; 
the standard deviation of the other wind and non-wind cases is 
similar. The maximum shell expansion velocity is indicated by the 
grey solid line for 0.1 Myr. 

turbulence are not likely to explain the majority of ob¬ 
served features, most of which do show a progenitor with 
high confidence. A few of the observed shells also display 
evidence of interior ionized gas and warm dust, which re¬ 
quires a bright source. 

5. CONCLUSIONS 

We use magnetohydro dynamic simulations to inves¬ 
tigate the impact of stellar winds on molecular cloud 
structure and kinematics. We compute the CO( 1 - 0 ) 
cloud emission and compare it with observed data. We 
demonstrate that the morphology of the resulting sim¬ 
ulated shells is similar to structures observed in nearby 
clouds, which are attributed to spherical mass-loss from 
young stars. We explore a range of wind mass-loss rates 
and conclude that the observed size of the shells and 
their properties can only be created by progenitors with 
mass-loss rates of 10 “^ — lO“^M 0 yr“^, which is con- 
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sistent with the rates estimated from the observed CO 
molecular line data. This highlights a discrepancy of sev¬ 
eral orders of magnitude between observations and mod¬ 
els of wind launching, which predict mass-loss rates of 
< lO“^M 0 yr“^ for the stellar types identified. 

We evaluate the ability of the winds to drive super¬ 
sonic turbulence on scales of a few parsecs. A significant 
contribution to the turbulent cascade from winds hinges 
on the creation of non-local perturbations; otherwise the 
winds only deposit energy and momentum in the shells 
and do not offset global turbulent dissipation. While 
we find that the winds do create significant gas velocity 
enhancements compared to the initial cloud turbulence, 
which are especially notable in CO line emission, this 
appears to be a purely localized effect. Expansion slows 
considerably after the first 0.05Myr and most shells be¬ 
come pressure confined within a few O.lMyr. We find 
that the winds leave no signatures in the density and 
momentum distribution functions. The winds do create 
a feature at ^ Ipc in the velocity Fourier spectrum, but 
no corresponding feature appears in the power spectra of 
the integrated CO emission. We find no evidence, even in 
the strongest wind case, that the winds offset the global 
turbulent dissipation. 

Despite the localized nature of the wind deposition, we 
demonstrate that the energy and momentum of the shells 
can be considerable compared to the total cloud energy 
and momentum. However, energetics estimated using 
CO emission, which selects a particular range of interme¬ 
diate gas densities, tends to over-estimate the true wind 
impact. The high-velocity component {v > 1kms“^) of 
the CO emission suggests that 10-50% of the total cloud 
energy and momentum resides in shells. However, with 
the exception of the strongest wind case, the true contri¬ 
bution as measured using all the 3D gas information is 


not more than r\j 5%. 

We discuss alternative proposed origins for the shells, 
including turbulent voids, protostellar outflows, and 
wind variability, and we conclude that only variability 
presents a plausible alternative explanation. These re¬ 
sults have significant implications for our understanding 
of stellar winds, and our work suggests wind mass-loss 
rates are significantly underestimated by current mod¬ 
els. While the impact of winds from AB stars remains 
smaller than those of more massive O stars, AB stars are 
significantly more numerous and could provide signifi¬ 
cant impact in the aggregate. Consequently, underesti¬ 
mates of mass-loss rates have implications for molecular 
cloud properties and galaxy evolution within c osmolog- 
ical simulations (e.g., Agertz & Kravtsov||2014|). Future 
numerical work including gravity will examine the direct 
impact of wind feedback on star formation activity. 
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APPENDIX 

ANALYTIC COMPARISON OF WIND EXPANSION 

We test the wind implementation by computing the expansion of the wind launched by a single source embedded in 
a uniform density medium. We note the ORION hy drodynamic methodology has been benchmarked against a variety 
of standard tests i n pri or work (e.g., Li et al.||2Q12 ). 


Koo & McKee (1992), henceforth KM92, have solved for the expansion of wind bubbles in a uniform medium given 

in terms of the mechanical luminosity, L 


as 


some constant mass-loss rate. KM92 write the mass-loss rate, 

Min = 2L/v2 , where is the wind velocity. We rewrite and renormalize their solutions here in terms of M^. The 
size of an expanding radiative shell as a function of time is then (KM92, equation 3.1): 
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where nn = po/2.34 x 10“^^g is the Hydrogen number density. Figure 
of the wind shock front and the expected analytic solution. 
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shows good agreement between the location 


WIND EXPANSION IN A TURBULENT MACNETIZED MEDIUM 

Although the gas in molecular clouds is non-uniform due to turbulence, the KM92 solutions also apply to the average 
expansion with some modification. Since the winds here expand into a cold ISM, the swept-up shell material is also 
relatively cold (the entrained gas mass is much larger than the mass in the wind) and the winds can be described by 
the slow-wind (radiative) solution. 

Once the ram pressure of the shell is equal to the ambient pressure, the shell becomes “pressure-confined” and the 
expansion stops. The time that this occurs for a radiative bubble is (KM92, equation 3.8): 
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Fig. 12 .— Radius of the shell shock front as a function of time (blue diamonds) with ambient density po = 2.04 x 10 ^^gcm ^ (the 
fiducial cloud density in the mo dels) and a source with mass-loss rate ttIu) = 6 x 10“®MQyr“^. The grey solid line gives the corresponding 
analytic solution from equation |A1| 




Fig. 13. — Log density (left) and velocity (right) Fourier spectra computed using 256^ and 512^ grid resolutions for two different times. 
Power laws of k~^-^ and respectively, and approximate inertial range are denoted by grey lines. The spectra have been compensated 

by k?. 

where = Vwjcs is the effective Mach number of the wind as compared to the isothermal sound speed in the ambient 
medium. Given a wind velocity of 200 kms“^ and a gas temperature of 10 K, Cg — 0.2km s“^ and Mw = 1000. This 
predicts that a wind with IO^^Mq yr“^ would not become pressure confined for more than 3 Myr. However, the ambient 
medium here is both magnetized and turbulent. Consequently, it is appropriate to replace Cg with the effective sound 
speed, Cs,eff = + v \, where (Jid is the rms ID turbulent velocity dispersion and va is the Alfven speed. For 

our fiducial simulation, Ceff = 1.43kms“^, which means that a source with a mass-loss rate lO“^M 0 yr“^ will become 
pressure confined in tp = 0.068 Myr. 


FOURIER SPECTRA RESOLUTION COMPARISON 

The basegrid for the calculation that is used to genera te the initial turbule nt conditions is 256^. Larger grid sizes 
are generally preferred for computing turbulent statistics (Kritsuk et al. l2007t . Figure [^compares the density power 
spectrum for the 256^ data and for a fiattened grid of 512^ that takes into account the AMR level 1 data. In the latter 
case, strong shocks as well as the wind shells are captured at finer resolution. The power law exponents of the lower 
resolution spectra are ^ 0.1 fiatter than the higher resolution power laws, where the higher resolution is more similar to 
prior published values for turbulent calculations. The 256^ resolution momentum and velocity Fourier spectra power 
law fits are also fiatter by ~ 0 . 1 . 


Agertz, O., & Kravtsov, A. V. 2014, ArXiv e-prints, 
arXiv:1404.2613 
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